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ABSTRACT 


A method  for  numerically  generating  boundary- 
fitted  coordinate  systems  for  three-dimensional 
regions  containing  ship-like  bodies  is  the  subject 
of  this  report.  This  procedure  involves  a trans- 
formation which  maps  the  region  of  interest  in 
physical  space  onto  a region  in  computational 
space  where  a uniform  grid  is  defined.  The  result 
is  a curvilinear  coordinate  system  in  the  physical 
region  having  coordinate  lines  coincident  with  all 
boundary  contours.  The  finite-difference  solution 
of  a system  of  partial  differential  equations 
corresponding  to  a physical  problem  can  be  obtained 
on  the  regular  mesh  in  the  computational  range  space. 
Alternatively,  the  boundary-fitted  mesh  in  the 
physical  domain  space  may  be  used  as  the  basis  of  a 
finite-element  scheme  for  this  solution.  The 
mapping,  which  is  a generalization  of  conformal 
mapping,  is  found  as  the  solution  of  a system  of 
elliptic  partial  differential  equations  with  appro- 
priate boundary  conditions.  A general  description 
of  the  technique  is  presented  along  with  the 
mathematical  formulation  of  the  underlying  boundary- 
value  problem.  Also  included  are  methods  of 
transformation  control  which  are  needed  for  developing 
grid  systems  suitable  for  numerical  solutions  of 
problems  in  fields  such  as  fluid  dynamics.  Examples 
of  several  transformations  are  discussed  in  terms  of 
their  associated  physical  and  computational  regions 
and  the  elliptic  systems  used  for  their  generation. 


INTRODUCTION 

The  increasing  usage  of  boundary-fitted  (surface-oriented) 
coordinates  in  the  field  of  fluid  dynamics  is  indicative  of  this  tech- 
nique's value  as  a numerical  problem-solving  tool.  Typically,  the 
investigator  wishes  to  study  the  flow  of  a fluid  in  a particular  region 
by  analyzing  a partial  differential  equation  subject  to  certain  boundary 
conditions.  The  proper  choice  of  coordinate  system  is  often  crucial  to 
the  successful  solution  of  the  problem.  As  a result,  much  use  has  been 
made  of  "natural"  coordinate  systems  such  as  cylindrical  and  spherical 
coordinates,  although  the  practical  usefulness  of  such  systems  is 
restricted  to  a few  specialized  cases.  For  any  specific  geometry,  the 
method  of  boundary-fitted  coordinates  creates  a system  in  which  the 
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surfaces  themselves  are  coordinate  lines. 

The  procedure  involves  a numerically  generated  transformation  which 

maps  the  physical  region  under  consideration  onto  a simpler  computational 

region  where  a uniform  mesh  has  been  defined.  In  effect  two  grid  systems 

are  produced  by  this  transformation:  the  boundary-fitted  mesh  in  the 

domain  space  and  the  regular  mesh  in  the  range  space.  Either  of  the 

meshes  may  be  used  to  compute  the  numerical  solution  to  a partial 

differential  equation  corresponding  to  a physical  problem. 

In  obtaining  this  solution,  a finite-difference  approach  must  take 

into  account  variations  in  the  mesh  upon  which  the  computations  are  to  be 

done.  The  numerical  scheme  is  often  complicated  by  these  variations 

which  normally  occur  in  zones  where  increased  accuracy  is  desired  or 

where  irregularly-shaped  boundaries  are  resolved.  If,  however,  the 

calculations  are  performed  on  the  regular  mesh  in  the  computational  space, 

these  complications  disappear  since  the  mesh  spacing  is  uniform  there. 

This  usually  leads  to  a simpler  numerical  treatment  even  though  the 

equation  of  ultimate  interest  may  be  somewhat  altered  by  the  transfor- 
1* 

mat ion . 

Calculations  may  also  be  carried  out  in  the  physical  region  with 

the  grid  points  there  being  taken  as  the  nodes  of  a finite-element  mesh. 

The  important  aspect  here  is  not  the  transformation  itself,  but  the  fact 

that  the  network  of  mesh  points  in  the  physical  space  conforms  to  the 

boundary  contours.  Used  in  this  manner,  the  numerically  generated 

coordinate  technique  provides  an  efficient,  automatic  means  of  producing 
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a boundary-fitted  grid  system  suitable  for  finite-element  calculations. 

This  numerical  mapping  technique  was  extended  to  three-dimensional 

3 4 

domains  in  earlier  work  by  Ghia  and  Mastin.  The  purpose  of  the  present 
research  is  to  present  some  examples  of  new  types  of  transformations  that 
are  proving  useful  in  the  study  of  fluid  flow  about  ship-like  bodies. 

As  expected,  the  generalization  from  two  to  three  dimensions  consti- 
tutes a considerable  increase  in  the  complexity  of  the  entire  process. 

The  increased  computer  time  and  memory  requirements  have  been  overcome 
to  some  extent  through  the  use  of  the  Texas  Instruments  Advanced 

*A  complete  listing  of  references  is  given  on  page  17. 
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Scientific  Computer  (TIASC)  at  the  Naval  Research  Laboratory.  Since  the 
number  of  possible  mapping  configurations  also  rises  with  the  number  of 
dimensions,  the  experienced  user  of  this  approach  should  be  able  to 
produce  a coordinate  system  suitable  for  practically  any  flow  geometry  of 
interest . 

FINITE-DIFFERENCE  APPLICATION  TO  FLOW  PROBLEMS 

The  body-fitted  coordinate  technique  has  several  features  which  make 
it  attractive  for  use  in  conjunction  with  a finite-difference  scheme  for 
fluid  flow  problems.  One  of  these  features  is  the  three-dimensional 
mapping  capability  as  opposed  to  complex  variable  methods  which  are 
limited  to  two-dimensional  flow  problems.  Although  the  numerically 
generated  transformations  are  not  conformal  (the  orthogonality  of  coor- 
dinate lines  is  not  preserved),  there  is  no  real  problem  since  orthogon- 
ality is  not  required  at  boundaries  to  obtain  an  accurate  representation 

S 

of  the  normal  derivative. 

Another  important  property  is  that  all  surface  contours  in  the 
physical  region  are  coincident  with  coordinate  surfaces.  This  allows  for 
the  accurate  representation  of  boundary  conditions  regardless  of  the  shape 
and  number  of  bodies  present. 

All  calculations,  both  to  generate  the  coordinate  system  and  to  solve 
the  equations  governing  the  flow  problem,  are  done  on  the  fixed,  uniform 
mesh  in  the  transformed  region.  Since  the  mesh  must  be  recomputed  when- 
ever the  boundaries  change  or  deform,  the  coefficients  of  the  difference 
operators  may  also  change.  However,  this  updating  of  the  transformation 
required  a relatively  small  portion  of  the  total  computing  time  in  a 
recent  study  of  two-dimensional  unsteady  free  surface  flow.^ 

Also,  coordinate  lines  in  the  physical  region  may  be  concentrated  in 
areas  where  more  resolution  or  higher  accuracy  is  desired.  This  concen- 
tration or  control  of  the  coordinate  lines  can  be  accomplished  in  several 
ways:  1)  by  changing  the  mapping  configuration  itself,  2)  by  varying  the 
underlying  elliptic  generating  equations,  and  3)  by  altering  the 
distribution  of  grid  points  on  the  physical  boundaries.  A discussion  of 
coordinate  system  control  is  found  in  a later  section. 
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MATHEMATICAL  FORMULATION 

The  technique  of  numerical  coordinate  system  generation  for  two- 
dimensional  regions  and  a finite  difference  scheme  for  its  implementation 
were  presented  in  a previous  publication.^  Here  we  develop  the  three- 
dimensional  counterpart  in  a slightly  different  manner  using  tensor 
notation. 

A three-dimensional  region  of  arbitrary  shape,  R,  is  to  be  trans- 
formed into  the  rectangular  region,  R' , as  shown  in  Figure  1.  For 
convenience,  let  xi  be  the  usual  Cartesial  coordinates  and  ui  be  the 
transformed  coordinates.  The  u*  are  obtained  as  solutions  of  the  system 

V^1  = Pt,  i = 1,2,3  (1) 

where  V2  is  the  Laplacian  operator  in  Cartesian  coordinates,  so  that 


72  = 


Sx'Sx1  8x23x2  3x33x3 
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where  S , S„,  ...,  S are  the  surfaces  indicated  in  Figure  1. 
1 Z D 
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There  are  several  aspects  of  this  general  boundary-value  problem 
which  are  somewhat  arbitrary  and  must  be  specified  by  the  user.  These 
are  the  shape  of  the  computat.  lal  region,  the  in  Equation  (2)  and 
the  u^  in  Equations  (3).  This  specification  can  be  made  in  a manner  that 
will  greatly  influence  the  properties  of  the  resulting  coordinate  system. 
The  ability  to  control  the  coordinate  system  which  is  generated  is  an 
important  topic  and  is  the  subject  of  the  next  section. 

Since  all  computations  are  to  be  done  in  the  transformed  region  R' , 
we  must  interchange  the  roles  of  the  dependent  and  independent  variables 
in  Equation  (1) . The  inverted  system  is 


Xp  S^x1 
8 3uX3u^ 


+ P ~o, 

v 3uv 


i = 1,2,3 


(4) 


where  g^u  are  the  components  of  the  contravar iant  tensor  associated  with 
the  Euclidean  metric  tensor  g . The  transformed  boundary  conditions  are 

Ml 

simply  the  physical  coordinates  of  the  mesh  points  on  S,  , S , ...  , S . 

i Z o 

Equation  (4)  is  approximated  using  second-order,  central 
differences  for  all  derivatives  involved,  and  the  resulting  difference 
equations  solved  on  the  uniform  mesh  in  the  computational  region.  The 
difference  equations  and  iteration  scheme  are  discussed  in  detail  for  the 
two-dimensional  system  in  earlier  publications^ ’^Accelerated  Gauss-Seidel 
iteration  was  used  to  produce  the  examples  in  this  report,  although  this 
may  not  be  the  most  efficient  method  for  computers  with  vector  processing 
capability.  A vectorizable  iteration  method  which  sweeps  the  grid  points 
in  an  alternating  manner,  the  so-called  "red-black"  method,  was  used  by 
Haussling^  for  a two-dimensional  problem. 


COORDINATE  SYSTEM  CONTROL 

As  mentioned  earlier,  there  are  several  methods  of  coordinate 
system  control.  These  methods  can  be  used  by  the  researcher  to  tailor 
the  coordinate  system  to  his  particular  problem.  Although  there  is  a 
great  deal  of  flexibility  inherent  in  this  method,  care  must  be  exercised 
in  order  to  obtain  a proper  mapping. 
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The  first  way  of  influencing  the  coordinate  system  is  by  changing 
the  mapping  configuration  itself.  For  simplicity,  the  physical  domain  in 
Figure  1 was  taken  as  a deformed  cube.  If  R had  been  less  regular  in 
shape,  a computational  region  made  up  of  more  than  one  cube  might  have 
been  necessary  to  generate  a desirable  coordinate  system. 

The  second  area  of  flexibility  lies  in  the  elliptic  equations  that 
are  taken  as  the  basis  of  the  generating  system.  Although  any  equation 
exhibiting  an  extremum  principle  may  be  used,  the  Poisson  equation  seems 
to  suffice  in  most  cases  since  we  have  the  in  Equation  (2)  at  our 
discretion.  The  ability  to  vary  these  non-homogeneous  terms  provides  us 
with  an  effective  means  of  controling  the  distribution  of  the  coordinate 

lines  in  the  physical  region  R. 

17  8 

In  several  earlier  works,  ’ ' the  inhomogeneous  terms  were  chosen 
such  that  the  coordinate  lines  were  attracted  to  or  repelled  from  certain 
specified  points  in  the  physical  region.  This  approach,  while  eifective, 
proved  to  be  both  tedious  and  time-consuming.  Often,  the  procedure 
consisted  of  making  an  initial  estimate  of  the  parameters  in  the  P. 
terms,  generating  a mesh  with  these  parameters,  and  displaying  it 
graphically.  This  process  was  repeated,  with  slight  changes  in  the 

control  parameters,  until  a suitable  transformation  was  produced. 

9 

Thompson  presented  a method  of  choosing  the  inhomogeneous  terms  so 
that  the  spacing  of  the  coordinate  lines  on  the  boundary  is  maintained 
throughout  the  entire  region.  This  choice  of  source  terms  was  effective 
in  creating  meshes  for  problems  involving  free  surfaces  with  small 
elevations.  Large  amplitude  waves,  however,  are  not  resolved  satisfac- 
torily by  this  Poisson  system.  A modification  of  Thompson's  derivation 
leads  to  the  following  definition  of  the  P^  of  Equation  (2): 


P. 

l 


-32x1/3u13u 

4 3 

(3x1/3u1) 


c . 

l 


(5) 


where  u^  = c^  are  the  boundaries  on  which  the  x*  spacing  is  to  be 
preserved.  These  source  terms  have  produced  improved  coordinate  systems 
for  flow  regions  with  both  largo  anl  small  free  surface  elevations. 


The  third  means  of  coordinate  system  control  is  the  specification  of 
the  two  non-constant  coordinates,  in  Equations  (3),  on  the  boundaries. 
The  user  does  not  have  complete  freedom  here  since  this  input  characterizes 
the  boundaries,  but  in  most  cases  there  is  sufficient  leeway  to  provide  a 
coordinate  system  with  the  desired  properties  on  the  boundaries. 

It  is  well  to  note  that  often  the  most  satisfactory  coordinate  system 
is  obtained  when  these  three  methods  of  control  are  used  in  the  proper 
combination.  One  might  first  choose  a mapping  configuration  that  has  the 
desired  general  properties.  Some  properties  considered  might  be,  for 
example,  which  coordinate  will  be  constant  along  each  boundary,  how  many 
coordinate  lines  will  intersect  each  boundary,  and  the  shape  of  the 
computational  region,  etc.  After  the  general  configuration  has  been 
determined,  the  generating  equations  must  be  decided  upon.  A good  initial 
choice  seems  to  be  a system  of  Poisson  equations  with  source  terms  that 
preserve  the  boundary  spacing  throughout  the  region  as  described  above. 
Finally  one  can  adjust  the  distribution  of  grid  points  on  the  boundaries 
until  a suitable  mapping  is  achieved. 

EXAMPLES 

As  indicated  earlier,  the  flexibility  of  the  body-fitted  coordinate 
technique  is  one  of  its  most  important  assets.  There  is  no  restriction 
to  simple  regions,  although  multi-connected  and  other  complex  geometries 
may  require  a combination  of  simple  mappings.  The  examples  presented  in 
this  section  were  calculated  to  demonstrate  a few  of  the  many  possible 
transformations  that  can  be  produced. 

First,  we  consider  a physical  domain  which  contains  a thin  body 
intersecting  one  boundary.  This  example  may  be  thought  of  as  a model  of 
a ship  hull  in  a free  surface.  The  natural  transformation  for  a body  of 
this  shape  is  one  which  maps  the  object  to  a portion  of  a coordinate  plane. 
Figure  2 shows  the  body-fitted  coordinate  system  in  physical  space 
generated  with  zero  source  terms,  P^,  and  the  corresponding  uniform  mesh 
in  computational  space.  Two  cross-sectional  views  of  this  curvilinear 
coordinate  system  are  shown  in  Figure  3.  (Table  1 gives  a summary  of  the 
computer  requirements  for  Figures  2-6.) 

8 


1 


9 


In  the  next  example  we  generate  a grid  system  about  a circular 
cylinder  of  finite  length.  This  mapping  provides  us  with  a coordinate 
system  in  physical  space  which  is  quite  different  from  the  usual  cylindri- 
cal coordinates.  For  illustrative  purposes.  Figure  4 shows  only  the 
lower  half  of  the  physical  and  transformed  regions.  Again,  the  source 
terms,  P^,  are  zero. 

To  demonstrate  the  multi-body  capability,  we  consider  another  ship- 
like form  in  conjunction  with  two  smaller  objects  close  by.  In  generating 
this  single  transformation  we  must  take  into  account  three  separate 
bodies.  The  largest  of  the  three  has  a circular  cross-section  with  a 
rounded,  sloping  bow  and  a flat  stern.  This  portion  is  handled  in  a 
manner  similar  to  that  used  in  the  second  example.  The  two  smaller 
bodies  at  the  ship's  stern  are  thin  and  are  transformed  as  in  the  first 
example.  The  mesh  systems  in  physical  and  computational  space  generated 
with  zero  source  terms  are  shown  in  Figure  5. 

As  an  example  of  coordinate  line  control,  consider  a physical  region 
simulating  a fluid  bounded  by  straight  walls,  a bottom,  and  a free 
surface.  Figure  6 shows  two  coordinate  systems  generated  for  such  a 
physical  region:  one  generated  with  zero  source  terms  and  one  generated 
with  source  terms  as  given  in  Equation  (5) . We  see  that  the  non- 
homogeneous  system  has  produced  a mesh  with  a more  desirable  configuration 
near  the  curved  surface. 

It  is  appropriate  at  this  point  to  mention  the  role  of  computer 
aided  graphics  in  numerical  mesh  generation.  Usually,  the  process  of 
obtaining  a suitable  coordinate  transformation  for  a specific  problem  is 
one  of  trial  and  error.  Interactive  plotting  routines  are  useful  in 
shortening  the  time  needed  by  helping  to  locate  errors  in  input  boundary 
data  and  by  displaying  the  final  output.  The  use  of  these  automated 
graphing  techniques  has  proven  to  be  practically  indispensible  in  the 
analysis  of  numerically  generated  coordinate  systems.  This  is  especially 
true  for  three-dimensional  systems.  The  figures  in  this  report  were 
generated  by  IMAGE, ^ an  interactive  data  display  package  available  in 
the  Computation,  Mathematics,  and  Logistics  Department  of  DTNSRDC.  For 
details  on  the  use  of  this  and  other  similar  programs,  contact  Code  1843. 
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Figure  5 - Multi-Body  Transformation 
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TABLE  1 - COMPUTER  REQUIREMENTS  FOR  FIGURES 


Figures  2 
and  3 

Figure  4 

"1 

Figure  5 

Figure  6 

Total  Number 
of  Points 

2541 

6629 

7553 

9261 

Memory 

Requirements 

131  K. 

204  K 

222  K 

236  K 

Computer 

Used 

TIASC 

CDC 

6400 

CDC 

6400 

CDC 

6400 

CPU 

Time 


7.1  sec 


20.9  sec 


22.3  sec 


34.2  sec 


FUTURE  PLANS 


Three-dimensional  versions  of  the  boundary-fitted  coordinate 
generating  programs  are  operational  on  the  CDC  6000  series  computers  at 
DTNSRDC  and  on  the  TIASC  at  the  Naval  Research  Laboratory.  Work  is 
continuing  on  refining  these  programs  in  preparation  for  the  numerical 
solution  of  free  surface  and  ship  wave  problems. 

Problems  currently  under  investigation  include  1)  the  non-linear 
effects  of  water  waves,  2)  large  amplitude  waves  nearing  the  point  of 
breaking,  and  3)  ships  undergoing  a slamming  motion  in  a free  surface. 
Most  likely  these  and  other  large  fluid  flow  problems  will  be  solved  on 
the  TIASC  because  of  its  large  memory  and  vectorization  capability. 
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